This invention relates to a large-scale multiplication with addition operation (mul-add operation) method and system and in particular to a large-scale multiplication with addition operation method and system for efficiently performing recursive multiplication with addition operation.
In image processing and simulation of physical phenomenon, continuous integral or differential calculation is often replaced with a discrete form, for example, a recursion formula to obtain an approximate numeric solution. An example is a polynomial containing a plurality of terms different in degree expressed in the following expression (1) in actual calculation performed by hardware of a computer, etc.,: EQU y=a.sub.0 X.sup.n +a.sub.1 X.sup.n-1 +a.sub.2 X.sup.n-2 +. . . +a.sub.n-1 X+a.sub.n (1)
As a technique of efficiently performing the calculation of the expression (1), the expression (1) is represented in the recursive multiplication with addition form for replacement with simple multiply with add processing repetition. Recursive representation equivalent to the expression (1) is shown in expression (2): ##EQU1##
Thus, the value of the polynomial of expression (1) can be found efficiently by simple repetition of recursive multiplication with addition operation. Therefore, if hardware capable of performing simple multiplication with addition operation at high speed is used, the above-mentioned processing can be performed efficiently. In expression (2) shown above, multiplication with addition operation may be performed n times recursively to obtain the value of sn. In the description to follow, the number of times recursive multiplication with addition operation is performed corresponding to the degree of a polynomial will be referred to as recursive order n. In the examples in expressions (1) and (2), the recursive order matches the number of times recursive multiplication with addition operation is performed.
Some simulation in a chemical field requires recursive calculation in complicated form. An example is ab initio molecular orbital calculation using a variational method. The purpose of this ab initio molecular orbital calculation is to describe electron behavior in a molecule in first principle calculation without using any empirical parameters.
The ab initio molecular orbital calculation method is positioned as the foundation for molecule simulation and is used to analyze material structures and chemical reaction in detail as an industrially important technique. In the ab initio molecular orbital calculation method, the reciprocal of an exponential function having as an exponent, the distance between the atomic nucleus of each atom and each orbital electron making up a molecule multiplied by an empirical constant (which is called orbital exponent) is used as an atomic orbital, which will be hereinafter referred to as a basis function (or simply z basis), and a plurality of the basis functions are provided for one atom. The basis functions are linearly combined, thereby describing the wave function of the orbital electron in the molecule, namely, the molecular orbital. Determining coefficients for linear combination of the basis functions in the molecular orbital is the main processing content of the ab initio molecular orbital calculation method.
More particularly, spatial (electron spin is not considered) orbit .phi.i of electron in molecule can be approximated by linear combination of basis functions .chi. as shown in the following expression (3) and (2n electron system) molecular state (wave function) .psi. satisfying Pauli exclusion principle considering electron spin can be expressed by Slater determinant shown in the following expression (4): EQU .phi..sub.i =.SIGMA..sub..mu. .chi..sub..mu. C.sub..mu.i (3)
where
i: Number representing ith molecular orbital .phi. (i, j, k, l, . . . , each being number of molecular orbital .phi.)
.mu.: Number representing .mu.th basis function .chi. (.mu., .nu., .lambda., .sigma., . . . , each being number of basis function .chi.
.SIGMA..mu.: Sum total regarding .mu. (sum of N basis functions when the number of basis functions is N)
C.mu.i: Constant for linear combination ##EQU2## where
.alpha.(x): Upward state of xth electron spin, and
.beta.(x): Downward state of xth electron spin
Expected value of molecular energy e expressed in the following expression (5) can be found from the expression (4) and Hamiltonian H (=H1+H2, H1: one electron portion, H2: two electron portion) for 2n electron system: EQU .epsilon.=.intg..PSI.H.PSI.d.tau.=.intg.(H.sub.1 +H.sub.2).PSI.d.tau.=2.SIGMA..sub.i H.sub.i +.SIGMA..sub.i .SIGMA..sub.j (2J.sub.ij -K.sub.ij) (5)
where .SIGMA.i and E.SIGMA.j denote the sum total regarding n molecular orbitals, Hi denotes core integral, Jij denotes Coulomb integral, and Kij denotes exchange integral.
The expression (5) can be expressed by the following expression (6) at the basis function level and the integral expressed by [.mu..nu., .lambda..sigma.] is called electron repulsion integral (ERI): ##EQU3##
C.mu.i: Constant for linear combination, i=1-n, .mu.=1-N, n=number of molecular orbitals, N=number of atomic orbitals (bases)
Since the expected value of molecular energy .epsilon. in the expression (6) contains unknown linear combination coefficients C, a numeric value cannot be obtained. Hartree-Fock-Roothaan variational method (HFR method) shown in the following expression (7) is known as a method for defining a coefficient matrix having the linear combination coefficients C as elements so as to minimize the expected value .epsilon. and finding wave function .psi. in a ground state (stable state with minimum energy): ##EQU4##
where
F.mu..nu.: Fock matrix element, P.lambda..sigma.: Density matrix element
The Fock matrix elements and the density matrix elements have numeric values for 1 to n basis functions .mu., .nu., .lambda., .sigma. and are represented in the N.times.N matrix form. The coefficient matrix C can be found by solving the expression (7). Since the expression (7) exists for each of i=1-n and .mu.=1-N, n.times.N simultaneous equations result.
The density matrix P is used to calculate the coefficient matrix C in the above-mentioned expression. Since the density matrix P is calculated from the coefficient matrix C, specifically the density matrix P calculated by previously giving initial value to the coefficient matrix C is used to find a Fock matrix, then a new coefficient matrix C is found. Next, the process is repeated until the difference between the coefficient matrix C on which the density matrix P is based and the coefficient matrix C obtained as the calculation result becomes sufficiently small (self-consistent). If the coefficient matrix C is found, linear combination coefficient is found and molecular orbital can be found.
By the way, the calculation of the Fock matrix from the density matrix P requires the calculation amount and the storage capacity proportional to the fourth power of the number of basis functions. Thus, the ab initio molecular orbital calculation method is applied only to a small-scale molecule system of about 100 atoms under present circumstances. Development of a calculation system dedicated to ab initio molecular orbital calculation with application to a molecule system of more than 1000 atoms in view is indispensable to make molecular theoretical elucidation of life and chemical phenomena more realistic.
That is, a substantially large computation load is calculation of Fock matrix elements containing electron repulsion integral, because calculation for all .mu. and .nu. is essential and the sum regarding .lambda. and .sigma. must be calculated for combinations of .mu. and .nu.. Therefore, calculation of N.sup.4 electron repulsion integrals is required per self-consistent iterative calculation.
An algorithm suitable for high-speed calculation for the ab initio molecular orbital calculation is shown in document A (S. Obara and A. Saika, "Efficient recursive computation of molecular integrals over Cartesian Gaussian functions," Journal of Chemical Physics, vol.84, no.7, pp.3963-3974, April 1986). In this algorithm, basis function .chi. is expressed by the following expression (8): EQU .chi.=.chi.(r, .zeta., n, R)=(x-R.sub.x).sup.nx (y-R.sub.y).sup.ny (z-R.sub.z).sup.nz e x p {-(r-R).sup.2 } (8)
where r, n, and R are vectors, r is electron coordinates having components (x, y, z), .zeta. is an orbital exponent (spatial extent of wave function), n is an angular momentum index having components (nx, ny, nz), and R is atomic nucleus coordinates at the center of atomic orbital having components (Rx, Ry, Rz). The sum total of the angular momentum index n (nx+ny+nz) represents the angular momentum; when the sum total is 0, it is called s orbit, when 1, p orbit, and when 2, d orbit, . . . To calculate a molecular orbital, a plurality of combinations (n, .zeta.) are provided for each of atoms constituting the molecule whose orbital is to be calculated and all sets are used as a basis function set. One basis function can be expressed by combination (n, .zeta., R).
As described above, the basis function .mu., .nu., .lambda., .sigma. is expressed by the angular momentum index, the orbital exponent, and the central atomic nucleus coordinates. However, in the description that follows, for simplifying the representation, angular momentum indexes a, b, c, and d will be used to express electron repulsion integral as [ab, cd] in place of the basis functions .mu., .nu., .lambda., and .sigma..
In the algorithm, electron repulsion integral occupying most of calculation for determining linear combination coefficient is expressed by a recursion formula in the form of the sum of terms consisting of products of auxiliary integral and coefficients, shown in the following expression (9) (see expression (39) in the document A): ##EQU5##
where i: Any of x, y, and z components
Similar recursion formulae also hold for [a (b+1i), cd], [ab, (c+1i)d], and [ab, c (d+1i)] and resemble expression (9), but will not be described. In the description to follow, expression (9) and its similar expressions may be collectively called expression (9).
In the above-shown expression (9), [.cndot..cndot..cndot.].sup.(m) represents auxiliary integral, wherein a, b, c, and d are represented as S, p, d, f, . . . corresponding to angular momentum 0, 1, 2, 3, . . . and (m) is an auxiliary index. Auxiliary integral with m=0 matches electron repulsion integral. Further, (a+1i) means that i component of angular momentum a (i is one of x, y, and z) is incremented by one; (a-1i) means that i component of angular momentum a is decremented by one. In addition, Pi, Ai, Wi, .zeta., p, and .eta. are constants and Ni (a) and the like are variables dependent on a and i, which are similar to those in the document A and therefore will not be discussed in detail here.
Paying attention to the angular momentum of the leftmost basis function of the four basis functions constituting auxiliary integral, at least one is decreased in eight terms on the right-hand side than the left-hand side by recursion formula expansion using expression (9). If the angular momenta of the four basis functions constituting auxiliary integral are all 0, namely, for [ss, ss].sup.(m), expression (10) shown below (see expression (44) in the document A) can be used to obtain a numeric value. Therefore, electron repulsion integral constituted of the four basis functions having angular momenta of any magnitude can be expressed using auxiliary integral consisting of only basis functions each finally having the angular momentum 0 after expression (9) is repeatedly applied for decreasing the angular momentum; the numeric value of electron repulsion integral can be easily obtained. EQU [0.sub.a 0.sub.b,0.sub.c 0.sub.d ].sup.(m) =(.zeta.+.eta.).sup.1/2 K(.zeta..sub.a,.zeta..sub.b, A, B)K(.zeta..sub.c,.zeta..sub.d,C, D)F.sub.m (T) (10)
Specifically, electron repulsion integral is expressed in the form of the sum of products of eight or less auxiliary integrals and coefficients according to expression (9) and the angular momentum is decreased further according to expression (9) for the auxiliary integrals. This procedure is repeated and all the procedure to auxiliary integral with all angular momenta 0 is recorded. Next, the recorded procedure is reversed starting at the auxiliary integral with all angular momenta 0 having a known numeric value for finding the numeric value of each auxiliary integral and finally obtaining the numeric value of electron repulsion integral.
As an example, the electron repulsion integral [pxpx, pxpx] calculation procedure can be expanded into the following seven recursion formulae of expressions (11)-(17) from expression (9), etc.,: ##EQU6##
The electron repulsion integral [pxpx, pxpx] is [pxpx, pxpx].sup.(0) in auxiliary integral expression and can be expressed in the form of the sum of products of six auxiliary integrals [pxpx, pxs].sup.(0), [pxpx, pxs].sup.(1), [spx, pxs].sup.(1), [pxs, pxs].sup.(1), [pxpx, ss].sup.(0), and [pxpx, ss].sup.(1) and coefficients using expression (11). The auxiliary integrals [pxpx, pxs].sup.(0) and [pxpx, pxs].sup.(1) can be expanded into recursion formulae using expression (12), the auxiliary integral [spx, pxs].sup.(1) can be expanded into a recursion formula using expression (13), the auxiliary integral [pxs, pxs](.sup.1) can be expanded into a recursion formula using expression (14), and the auxiliary integrals [pxpx, ss].sup.(0) and [pxpx, ss].sup.(1) can be expanded into recursion formulae using expression (15).
When a further description is given taking the auxiliary integral [pxpx, pxs].sup.(0) as an example, it can be expressed in the form of the sum of products of four auxiliary integrals [pxpx, ss].sup.(0), [pxpx, ss].sup.(1), [spx, ss].sup.(1), and [pxs, ss].sup.(1) and coefficients using expression (12). Using thus recursion expression expansion repeatedly, [pxpx, pxpx] can be expressed in the form of the product sum of [ss, ss].sup.(m) (where m=0-4) and coefficients finally using expression (16) or (17).
The numeric calculation procedure is executed in the reverse order to the recursion expansion. First, the numeric values of the auxiliary integrals [ss, ss].sup.(m) are provided for m=0-4 and the numeric values of the auxiliary integrals [pxs, ss].sup.(m) and [spx, ss].sup.(m) are found for m=0-3 using expressions (16) and (17). Similar operation is performed to find the numeric values of auxiliary integrals [pxpx, pxs].sup.(0), [pxpx, pxs].sup.(1), [spx, pxs].sup.(1), [spx, pxs].sup.(1), [pxs, pxs].sup.(1), [pxpx, ss].sup.(0), and [pxpx, ss].sup.(1) are found and finally expression (11) can be used to obtain the numeric value of the electron repulsion integral [pxpx, pxpx].
Relating the final numeric value y=sn in expression (2) to the auxiliary integral [pxpx, pxpx].sup.(0) and s1 in expression (2) to the auxiliary integral [pxs, ss].sup.(O), for example, the above-described operation to find the numeric values of the auxiliary integrals [pxpx, pxpx].sup.(O) and [pxs, ss].sup.(0) is recursive calculation and the recursive order of the recursive calculation to find the numeric value of the auxiliary integral [pxpx, pxpx].sup.(0) is larger than the recursive calculation to find the numeric value of the auxiliary integral [pxs, ss].sup.(0).
The ab initio molecular orbital calculation algorithm is represented by a more complicated form as compared with the recursive calculation represented in expression (2) and the recursive order cannot be expressed simply by the number of times multiplication with addition operation is performed. In the description to follow, the numeric value of the maximum one of the angular momenta contained in four basis functions making up electron repulsion integral or auxiliary integral will be called recursive order in the ab initio molecular orbital calculation algorithm. For example, the recursive order to find the numeric value of [ss, ss] is 0 and the recursive orders to find the numeric values of [pxs, ss], [pxs, pxs], and [pxpx, pxpx] are all 1.
A special purpose system for calculation containing electron repulsion integral requiring an enormous amount of calculation wherein recursive calculation is enabled in ab initio molecular orbital calculation is proposed by Shirakawa et al. (Choukousoku bunshikidou keisansenyouki MOEno architecture, Denshi Jyouhou Tuushingakkai gijyutu houkoku, CPSY96-46 (1996-05)). This multiplication with addition operation system comprises a host computer 1, a bus 2, a large-capacity storage unit 3, and processor units 4-0, 4-1, 4-2, . . . , as shown in FIG. 10. The host computer 1, the large-capacity storage unit 3, and the processor units 4-0, 4-1, 4-2, . . . can communicate with each other through the bus 2.
Calculation in the multiplication with addition operation system will be discussed with reference to FIG. 6. First, the host computer 1 prepares base information based on specification by the user at step 200. The base information consists of parameters such as the position coordinates of the atomic nucleus of each of the atoms constituting a molecule, types of atoms, orbital exponents describing basis functions of an orbital electron in each type of atom, angular momentum, and the like. In the description, the number of basis functions is N. The host computer 1 transmits the base information via the bus 2 to the processor units 4-0, 4-1, 4-2, . . . at step 202. It sets a trial value (initial value) of the linear combination coefficient of the basis functions at step 204 and transmits density matrix elements Pij (i=1-N and j=1-N) calculated from the trial value to the large-capacity storage unit 3 at step 206. Next, the host computer 1 determines basis function index assignment for each processor unit at step 208 and transmits it to each processor unit at step 210. The basis function index is provided for specifying the first one of the four basis functions constituting electron repulsion integral.
In the description to follow, electron repulsion integral [rs, tu] will be expressed simply as grstu (r, s, t, u=1-N). Calculation of N.sup.3 electron repulsion integrals gRstu (s, t, u=1-N) corresponding to one r=R according to the basis function index, is assigned to one processor unit.
The processor unit calculates electron repulsion integrals gRstu (s, t, u=1-N) for all s, t, and u for assigned r=R and calculates a part of Fock matrix element, F'Rs (s=1-N), according to the following expression (18): EQU F'.sub.rs =.SIGMA.P.sub.tu (g.sub.rstu -g.sub.rtsu /2) (18)
In expression (18), .SIGMA. means the sum total regarding all t and u and Ptu is a density matrix element stored in the large-capacity storage unit 3. Before calculating one electron repulsion integral, the processor unit previously receives necessary density matrix elements from the large-capacity storage unit 3.
Upon completion of calculation of all electron repulsion integrals, the processor unit transmits F'Rs (s=1-N) to the host computer 1. Upon completion of calculation of all F'rs executed by the processor units according to reception from the host processors at step 212, the host computer 1 adds a core integral amount to the calculation results to find Fock matrix elements at step 214 and calculates linear combination coefficients of basis functions based on the Fock matrix made up of the Fock matrix elements at step 216. Since calculation of the linear combination coefficients uses density matrix elements calculated using linear combination coefficient predetermined as initial value, calculation is repeated until they become self-consistent as described above at step 218. When they become self-consistent, the linear combination coefficient is determined to be the calculation result at step 220. F'rs is only a part of Fock matrix element as described above, but hereinafter will be called Fock matrix element for convenience.
Next, the configuration of each processor unit and electron repulsion integral calculation in the processor unit in the system in FIG. 10 will be discussed. FIG. 11 shows the configuration of the processor unit. The processor unit in FIG. 11 comprises a bus interface circuit 11, a program memory 12, a control circuit 13, an instruction memory 14, a data memory 15, and a multiplication with addition operation circuit 16. The bus interface circuit 11 is connected to the bus 2 and enables communication between the processor unit and the host computer 1 or the large-capacity storage unit 3.
A program transmitted from the host computer 1 through the bus 2 is written into the program memory 12 via the bus interface circuit 11 only once at the molecular orbital calculation starting time. The control circuit 13 mainly generates a multiplication with addition operation procedure in accordance with the program in the program memory 12 and writes the procedure into the instruction memory 14. It also determines the electron repulsion integral calculation sequence, generates a coefficient used for product sum calculation, stores the coefficient in the data memory 15, controls the multiplication with addition operation circuit 16, and issues an address to the data memory 15. The multiplication with addition operation procedure generated by the control circuit 13 is written into the instruction memory 14 and base information transmitted by the host computer 1 to the processor unit, the coefficient generated by the control circuit 13, the auxiliary integral calculation result, and the like are written into the data memory 15. Density and Fock matrix elements are also written into the data memory 15 and their numeric values can be transferred between the data memory 15 and the host computer 1 or the large-capacity storage unit 3 via the bus interface circuit 11 as required. The multiplication with addition operation circuit 16 loads numeric values stored in the data memory 15, performs multiplication with addition operation, and stores the computation result in the data memory 15 according to the multiplication with addition operation procedure stored in the instruction memory 14 under the control of the control circuit 13.
Next, calculation of Fock matrix element F'Rs executed by the processor unit will be discussed with reference to FIG. 3. Upon reception of basis function index R from the host computer 1 at step 100, the control circuit 13 determines the calculation sequence of the electron repulsion integrals to be calculated at step 102. Next, it generates a multiplication with addition operation procedure to obtain numeric values according to the sequence and stores the procedure in the instruction memory 14 at step 104. The same multiplication with addition operation procedure can be used if the electron repulsion integrals whose numeric values are to be found have the same angular momentum combination. Basis function.chi. can be expressed by the following expression (19) like the above-shown expression (8): EQU .chi.=.chi.(L, .zeta., X) (19)
where X is position coordinates of the central atomic nucleus of basis function X, .zeta. is an orbital exponent indicating the spatial extent of the basis function .chi., and L is the angular momentum of the basis function .chi..
The recursion formula shown in expression (9) indicates that if angular momentum L combinations match, the numeric values of electron repulsion integrals can be obtained by executing the same multiplication with addition operation procedure with change only in the numeric values of coefficients with different position coordinates X and orbital exponents .zeta.. Therefore, to determine the electron repulsion integral calculation sequence, all electron repulsion integrals to be calculated are classified into types according to L combination and are calculated for each type.
Next, a method of calculating one of the classified electron repulsion integrals will be discussed. First, the control circuit 13 generates a product sum calculation procedure corresponding to the product sum procedure type and stores the product sum calculation procedure in the instruction memory 14 at step 104. Next, the control circuit 13 selects the electron repulsion integral grstu to be calculated from among the electron repulsion integrals in the product sum calculation procedure type at step 106.
This step 106 is equivalent to determining of s, t, and u because R is fixed. Next, the control circuit 13 loads density matrix elements Ptu and Psu from the large-capacity storage unit 3 into the data memory 15 at step 108. The density matrix elements Ptu and Psu are numeric values required to calculate Fock matrix element F'Rs according to expression (18) when the numeric value of gRstu is obtained. Thus, when calculation of the electron repulsion integral is started, the control circuit 13 may request the large-capacity storage unit 3 to send the data for reception of the data in the processor unit by the time the calculation of the electron repulsion integral is complete. The control circuit 13 generates coefficients and numeric values of auxiliary integrals [ss, SS].sup.(m) used for the multiplication with addition operation based on the basis function R, s, t, u information (.zeta., X) stored in the data memory 15 and stores them in the data memory 15 at step 110.
Next, the control circuit 13 outputs a multiplication with addition operation start instruction to the multiplication with addition operation circuit 16 at step 112. When receiving the start instruction, the multiplication with addition operation circuit 16 performs multiplication with addition operation while reading necessary data from the data memory 15 according to the multiplication with addition operation procedure written in the instruction memory 14, and stores the computation result in the data memory 15. If it is determined at step 114 that execution of the multiplication with addition operation procedure written in the instruction memory 14 is complete, the calculation of one electron repulsion integral gRstu is complete. Upon completion of the calculation of one electron repulsion integral gRstu, the control circuit 13 uses the multiplication with addition operation circuit 16 to update the values of F'Rs and F'Rt at step 116.
Next, if it is determined at step 118 that a electron repulsion integral not yet calculated is contained in one multiplication with addition operation procedure type, control returns to step 106 at which s, t, and u are changed for calculating the electron repulsion integral not yet calculated without updating the multiplication with addition operation procedure stored in the instruction memory 14. On the other hand, if it is determined at step 118 that all electron repulsion integrals contained in one multiplication with addition operation procedure type have been calculated, control goes to step 120 at which whether or not calculation for all s, t, and u is complete is determined. If it is not complete, control returns to step 104 at which another type of multiplication with addition operation procedure is generated. If it is determined at step 120 that the calculation for all s, t, and u is complete, the calculation of F'Rs (s=1-N) is now complete.
Mention is made here of the size of data written into the instruction memory 14. For description, the multiplication with addition operation performed by the multiplication with addition operation circuit 16 is expressed as the following expression (20): EQU A3=A2+C.times.A1 (20)
where A1 and A2 are the numeric values of auxiliary integrals obtained so far, C is a coefficient, and A3 is the numeric value of auxiliary integral provided by the multiplication with addition operation. For the multiplication with addition operation circuit 16 to perform the above-described operation, the storage locations of the numeric values in the data memory 15 must be written in the instruction memory 14. That is, to find the numeric value of A3, the storage locations (addresses) at which A1, A2, and C are previously stored in the data memory 15 must be stored in the instruction memory 14. Thus, the instruction memory 14 requires a capacity of about eight bytes for one multiplication with addition operation.
If the maximum angular momenta of basis functions in molecular orbital calculation, etc., is 1, the electron repulsion integral with the amount of data to be written into the instruction memory 14 reaching the maximum is [pp, pp] In this case, the order of the recursive calculation for calculating the electron repulsion integral is 1. At the time, the necessary number of times multiplication with addition operation is executed until the numeric value of the electron repulsion integral [pp, pp] is obtained since starting at auxiliary integral [ss, ss].sup.(m) is 623. Thus, the capacity of the instruction memory 14 required for the multiplication with addition operation is about 5K bytes and the entire processor unit can be formed of a one-chip integrated circuit using a known SRAM. The entire processor unit is thus formed only of a one-chip integrated circuit, whereby information transfer is speeded up and high-speed multiplication with addition operation processing is enabled.
However, the capacity of the instruction memory for storing the multiplication with addition operation procedure is not considered for multiplication with addition operation that can be executed at high speed and in a short time in the above-described system; the executable computation is limited by the capacity of the instruction memory and multiplication with addition operation of large degree or order cannot be executed.
For example, in electron repulsion integral calculation required for molecular orbital calculation, calculation using basis functions with a large angular momentum cannot be executed. That is, considering practical utility of complicated computation such as molecular orbital calculation, at least 2 is required as the maximum value of the angular momenta of the basis functions. In this case, the electron repulsion integral with the amount of data to be written into the instruction memory 14 reaching the maximum is [dd, dd]. The order of the recursive calculation for calculating the electron repulsion integral is 2. Starting at auxiliary integral [ss, ss].sub.(m) to obtain the numeric value of the electron repulsion integral [dd, dd], the necessary number of times multiplication with addition operation is executed is 33936 and the capacity required for the instruction memory 14 becomes about 270K-bytes.
However, to form a one-chip integrated circuit containing the instruction memory 14 together with the control circuit 13 and the multiplication with addition operation circuit 16, a one-chip integrated circuit having an enormous capacity must be designed; it is substantially impossible. Thus, the capacity may be distributed for providing a separate instruction memory outside the chip, but access to the separate instruction memory occurs and the access time to the separate instruction memory would become extremely larger than the access time to the minute instruction memory 14 on one chip, resulting in a factor for hindering high-speed molecular orbital calculation.
A linear recursive operation system is proposed to perform calculation represented in recursive expression at high speed (refer to the Unexamined Japanese Patent Application Publication No. Hei 4-141769). However, the linear recursive operation system is applicable only to operation in a specific format called linear recursive operation and cannot be applied to calculation represented in a complicated recursive format such as ab initio molecular orbital calculation.
An auxiliary integral calculation method of molecular integral is proposed to efficiently calculate electron repulsion integral (refer to the Unexamined Japanese Patent Application Publication No. Sho 63-262758). In the art, a method of performing only calculation of simple auxiliary integrals [ss, ss].sup.(m) at high speed, but calculation containing basis functions with a large angular momentum, namely, complicated recursive calculation of large order cannot be executed.